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ABSTRACT 

We present numerical hydrodynamical evolutions of rapidly rotating rela- 
tivistic stars, using an axisymmetric, nonlinear relativistic hydrodynamics 
code. We use four different high-resolution shock-capturing (HRSC) finite- 
difference schemes (based on approximate Riemann solvers) and compare 
their accuracy in preserving uniformly rotating stationary initial configura- 
tions in long-term evolutions. Among these four schemes, we find that the 
third-order PPM scheme is superior in maintaining the initial rotation law 
in long-term evolutions, especially near the surface of the star. It is fur- 
ther shown that HRSC schemes are suitable for the evolution of perturbed 
neutron stars and for the accurate identification (via Fourier transforms) of 
normal modes of oscillation. This is demonstrated for radial and quadrupo- 
lar pulsations in the nonrotating limit, where we find good agreement with 
frequencies obtained with a linear perturbation code. The code can be used 
for studying small-amplitude or nonlinear pulsations of differentially rotat- 
ing neutron stars, while our present results serve as testbed computations 
for three-dimensional general-relativistic evolution codes. 

Key words: Hydrodynamics — relativity — methods: numerical — stars: 
neutron — stars: oscillations — stars: rotation 



1 INTRODUCTION 

The hydrodynamical evolution of neutron stars in nu- 
merical relativity has been the focus of many research 
groups in recent years (see e.g. Font et al. 1999; Miller, 
Suen & Tobias 1999; Mathews, Marronetti & Wilson 
1998; Shibata 1999a, Nakamura & Oohara 1998; Baum- 
garte, Hughes & Shapiro 1999; Neilsen & Choptuik 
1999). So far, these studies have been limited to ini- 
tially non-rotating starsQ. However, the numerical in- 
vestigation of many interesting astrophysical applica- 
tions, such as the rotational evolution of proto-neutron 

* After this work was completed, we received a preprint by 
Shibata 1999b, in which evolutions of approximate solutions 
of rotating stars are shown. 



stars and merged neutron stars or the simulation of 
gravitational radiation from unstable pulsation modes, 
requires the ability of accurate long-term evolutions 
of rapidly rotating stars. We thus present here an ex- 
tensive investigation of the nonlinear hydrodynamical 
evolution of rotating neutron stars (first results of this 
study were presented in Stergioulas, Font & Kokkotas 
1999). As we focus on the implementation and prob- 
lems of the hydrodynamical evolution of initially sta- 
tionary configurations in stable equilibrium, we do not 
evolve the spacetime. This approximation allows us to 
evolve relativistic matter for a much longer time than 
present coupled spacetime plus hydrodynamical evolu- 
tion codes. When studying the evolution of perturbed 
rotating neutron stars, the approximation of a static 
spacetime still allows for qualitative conclusions to be 
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drawn, since pulsations of neutron stars are mainly a 
hydrodynamical process. 

The rotational evolution of neutron stars can be af- 
fected by several instabilities (see Stergioulas 1998 for 
a recent review). If hot proto-neutron stars are rapidly 
rotating, they can undergo a dynamical bar-mode in- 
stability (Houser, Centrella & Smith 1994). When the 
neutron star has cooled to about 10 K after its forma- 
tion, it can be subject to the Chandrasekhar-Friedman- 
Schutz instability (Chandrasekhar 1970; Friedman & 
Schutz 1978) and it becomes an important source of 
gravitational waves. It was recently found that the 
I = m = 2 r-mode has the shortest growth time of 
the instability (Andersson 1998; Friedman & Morsink 
1998) and it can transform a rapidly rotating newly- 
born neutron star to a Crab-like slowly-rotating pul- 
sar within about a year after its formation (Andersson, 
Kokkotas & Schutz 1999; Lindblom, Owen & Morsink 

1998) . In this model, there are several important ques- 
tions still to be answered: What is the maximum am- 
plitude that an unstable r-mode can reach (limited by 
nonlinear saturation)? Is there significant transfer of 
energy to other stable or unstable modes via nonlinear 
couplings? (Owen et al. 1998). Does the r-mode evo- 
lution in a uniformly rotating star lead to differential 
rotation? (Spruit 1999, Rezzolla et al. 1999, Friedman 

1999) . Such questions cannot be answered by computa- 
tions of normal modes of the linearized pulsation equa- 
tions, but require nonlinear effects to be taken into ac- 
count. In addition, the computation of the frequency 
of r-modes in a rapidly rotating relativistc star may be 
easier to achieve in nonlinear numerical evolutions than 
in a perturbative computation. It is therefore desirable 
to develop the capability of full nonlinear numerical 
evolutions of rotating stars in relativity. The present 
implementation of the hydrodynamical evolution is a 
first step in this direction. 

For our study we have developed an axisymmet- 
ric numerical code in spherical polar coordinates which 
uses high-resolution shock-capturing (HRSC) finite- 
difference schemes for the numerical integration of 
the general relativistic hydrodynamic equations (see 
Ibanez & Marti 1999 for a recent review of applica- 
tions of HRSC schemes in relativistic hydrodynamics). 
Such schemes have been succesfully used before, in 
the study of the numerical evolution and gravitational 
collapse of nonrotating neutron stars in 1-D (Romero 
et al. 1996). An alternative approach, based on pseu- 
dospectral methods, has been presented by Gourgoul- 
hon (1991). Using our code in 1-D, we can accurately 
identify radial normal modes of pulsation up to high 
harmonic order. In 2-D the code is suitable for the hy- 
drodynamical evolution of rotating stars and the study 
of axisymmetric modes of pulsation. 

We find that the evolution of rotating stars with 
HRSC schemes can be subject to a significant numer- 
ical difficulty at the surface of the star. This difficulty 



arises due to the fact that one of the evolved hydro- 
dynamical variables, which corresponds to the angular 
momentum, has a sharp discontinuity at the surface. 
The result of a time-dependent evolution is a secular 
drift of the evolved configuration from the initial ro- 
tation law, near the surface. We investigate this prob- 
lem, using two different classes of HRSC schemes: To- 
tal Variation Diminishing (TVD, Harten 1984) and Es- 
sentially Non-Oscillatory (ENO, Harten & Osher 1987; 
Harten et al. 1987). We use both second and third- 
order variations of these schemes and find that, for a 
given resolution, the third-order schemes are superior 
in dealing with this problem. In addition, we show that 
all schemes are suitable for identiying normal modes of 
pulsation in a hydrodynamical evolution, via Fourier 
transforms of the evolved variables. 

In a subsequent paper we will present frequencies 
of quasi-radial and axisymmetric pulsations of rapidly 
rotating neutron stars. We further plan to study the 
pulsations of differentially rotating, hot proto-neutron 
stars and determine possible nonlinear couplings of dif- 
ferent modes of pulsation. 

The paper is organized as follows: In Section 2 
we describe the initial equilibrium configurations which 
will be evolved with the hydrodynamical code. In Sec- 
tion 3, we present the explicit form of the equations of 
general relativistic hydrodynamics as implemented in 
our code. Section 4 is devoted to the description of the 
numerical implemetation of the evolution schemes. Fi- 
nally, Section 5 presents 1-D and 2-D tests of the code 
with the various implemented schemes. 



2 INITIAL CONFIGURATIONS 

Our initial models are exact numerical solutions of 
rapidly rotating relativistic stars, having uniform an- 
gular velocity O. The metric of the stationary and ax- 
isymmetric spacetime in quasi-isotropic coordinates is 

ds 2 = — e 2v dt 2 + B 2 e~ 2v r 2 sin 2 6{d<f> — cudt) 2 

+e 2a {dr 2 + r 2 d6 2 ), (1) 

where v, B, a and ui are metric functions (Butterworth 
& Ipser 1976). In the non-rotating limit the above met- 
ric reduces to the metric of spherical relativistic stars in 
isotropic coordinates. We use dimensionless quantities 
by setting c = G — M = 1 . 

We assume a perfect fluid, zero-temperature equa- 
tion of state (EOS), for which the energy density is a 
function of pressure only. The following relativistic gen- 
eralization of the Newtonian polytropic EOS is chosen: 

V = ^ +1/JV (2) 
e = po + Np, (3) 

where p is the pressure, e is the energy density, po is the 
rest-mass density, K is the polytropic constant and N 
is the polytropic exponent. This form of a relativistic 
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polytropic EOS is convenient to use, since it coincides 
with the ideal fluid EOS used in many relativistic hy- 
drodynamical codes, in the case of an adiabatic evolu- 
tion. 

The initial equilibrium models are computed us- 
ing a numerical code by Stergioulas & Friedman 
(1995) which follows the Komatsu, Eriguchi & Hatchisu 
(KEH, 1989) method (as modified in Cook, Shapiro 
& Teukolsky, CST, 1992) with some changes for im- 
proved accuracy. In the KEH method three of the 
four general-relativistic, partial-differential field equa- 
tions are converted to integral equations using appro- 
priate Green's functions. The boundary conditions at 
infinity are thus incorporated in the integral evalu- 
ation. We use a compactified coordinate (introduced 
by CST) that allows one to integrate over the whole 
spacetime. The code has been shown to be highly 
accurate in extensive comparisons to other numer- 
ical codes (Nozawa et al. 1998). A public domain 
version is available at the following URL address: 
http : //www . gravity .phys .uwm.edu/Code/rns . 

The initial model is supplemented by a uniform, 
nonrotating "atmosphere" of very low density (typi- 
cally 10 -6 or less times the central density of the star). 
This is necessary for computing non-singular solutions 
of the hydrodynamic equations everywhere in the com- 
putational domain. After each time-step, we reset the 
density, pressure and velocity in the atmosphere to 
their initial values. The influence of the atmosphere is 
thus minimized and restricted to the grid-cells through 
which the surface of the star passes. The radial velocity 
of the atmosphere is only set to zero after each time- 
step, if it is negative. This avoids unwanted accretion 
of the material in the atmosphere onto the star and at 
the same time allows the star to expand, during radial 
pulsations. 



3 GENERAL RELATIVISTIC 

HYDRODYNAMIC EQUATIONS 

The equations of general relativistic hydrodynamics are 
obtained from the local conservation laws of density 
current J M and stress-energy T M " 



V M J" = 0, 
V M T M " = 

with 



J" 

rpfJ,l> 



POM , 



pg 



(4) 
(5) 

(6) 
(7) 



for a general EOS of the form p — p(p, e). Greek (Latin) 
indices run from to 3 (1 to 3). This choice of the 
stress-energy tensor limits our study to perfect fluids. 

In the previous expressions V M is the covariant 
derivative associated with the four-dimensional met- 
ric tt M is the fluid 4- velocity and h is the specific 



enthalpy 



h=l + £+-t- 



P_ 

A) 



(8) 



with e being the specific internal energy, related to the 
energy density e by 



e 

Po 



(9) 



With an appropriate choice of matter fields the 
equations of relativistic hydrodynamics constitute a 
hyperbolic system and can be written in a flux con- 
servative form, as was first shown by Martf, Ibanez 
and Miralles (1991) for the one-dimensional case. The 
knowledge of the characteristic fields of the system 
allows the numerical integration to be performed by 
means of advanced HRSC schemes, using approximate 
Riemann solvers (Godunov-type methods). The multi- 
dimensional case was studied by Banyuls et al. (1997) 
within the framework of the 3+1 formulation. Fur- 
ther extensions of this work to account for dynami- 
cal spacetimes, described by the full set of Einstein's 
non-vacuum equations, can be found in Font et al. 
(1999). Fully covariant formulations of the hydrody- 
namic equations (i.e., not restricted to spacelike ap- 
proaches) and also adapted to Godunov-type meth- 
ods, have been recently presented by Papadopoulos and 
Font (1999). 

In the present work we use the hydrodynamic 
equations as formulated in Banyuls et al. (1997). Spe- 
cializing for the metric given by Eq. ([!]), the 3+1 quan- 
tities read 

(10) 

2 sin 2 19, (11) 

(12) 
(13) 

n 2 6, (14) 

where a is the lapse function (the tilde is used to avoid 
confusion with the metric potential a), /3^, is the az- 
imuthal shift and lij is the 3-metric induced on each 
spacelike slice. 

The (axisymmetric) hydrodynamic equations are 
written as a first-order flux conservative system of the 
form 



a 


V 

= e , 




= -ujB 2 e 




2a 




= e , 




2 2a 


lee 


= re, 




= BV 2 " 



du daf 

dt dr 



+ 



86 



(15) 



expressing the conservation of mass, momentum and 
energy, where u, f r , f e and s are, respectively, the state 
vector of evolved quantities, the radial and polar fluxes 
and the source terms. More precisely, they take the 
form 

U = (D, S r , S$, S<j,,T), (16) 

r = (Dv r ,S r v r +p,Sev r ,S c ,v r ,(r+p)v r ), (17) 
f = {Dv e ,S r v e ,Sev e +p,S 4 ,v e ,(r + p)v e ). (18) 
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The source terms can be decomposed in the following 

way 



at at 



dr 

with 7 = det 7y and 



(19) 



5 |~T' 



(20) 



with j = r, 6 1 , 0. Quantities T„ v stand for the four- 
dimensional Christoffel symbols. The definitions of the 
evolved quantities in terms of the "primitive" variables 
w = (p , Vj,e) are 



D = p W, 



t = p hW 2 -p-D, 
where W is the relativistic Lorentz factor 
W ■ 



1 



(21) 
(22) 
(23) 

(24) 



with v 2 — 7ijU l u J . The 3- velocity components are ob- 
tained from the spatial components of the 4-velocity in 
the following way 

V= w + ^- (25) 



4 NUMERICAL SCHEMES FOR THE 
HYDRODYNAMIC EQUATIONS 

We turn now to the description of the numerical im- 
plemention of the equations presented in the previ- 
ous section. As we are investigating the accuracy of 
the evolution of rotating neutron stars with different 
schemes, we include here a detailed technical descrip- 
tion of the different numerical methods and tools we 
have implemented in our code. A point we emphasize 
is the use of high-order polynomial cell-reconstruction 
procedures, to achieve a better representation of the 
angular-momentum discontinuity at the surface of the 
star. 



4.1 A high-resolution shock-capturing 
algorithm 

Our hydrodynamical code performs the numerical in- 
tegration of system ([li]) using Godunov-type (HRSC) 
methods. In a HRSC scheme, the knowledge of the 
characteristic fields (eigenvalues) of the equations, to- 
gether with the corresponding eigenvectors, allows for 
accurate integrations, by means of either exact or ap- 
proximate Riemann solvers. These solvers, which con- 
stitute the kernel of our numerical algorithm, compute, 
at every interface of the numerical grid, the solution of 



local Riemann problems (i.e., the simplest initial value 
problem with discontinuous initial data). Hence, HRSC 
schemes automatically guarantee that physical discon- 
tinuities appearing in the solution, e.g., shock waves, 
are treated consistently (the shock- capturing property). 
HRSC schemes are also known for giving stable and 
sharp discrete shock profiles and for having a high order 
of accuracy, typically second order or more, in smooth 
regions of the solution. 

In our code we perform the time update of sys- 
tem ((i]§ from t n to t n+1 according to the following 
conservative algorithm: 



n+l n 
U i,j = U i,j 



At 
Ar 
At 



(fi+l/2,j — fj-l/2,j) 



~~ ^(S*,3 + l/2 _ gij-1/2) 
+ At Si j , 



(26) 



improved with the use of consecutive sub-steps to gain 
accuracy in time (see Shu & Osher 1989) . Index n rep- 
resents the time level and the time (space) discretiza- 
tion interval is indicated by At (Ar, A8). The "hat" in 
the fluxes denotes the so-called numerical fluxes which 
are computed by means of an approximate Riemann 
solver according to a generic flux-formula (supressing 
index j): 



f, 



+1/2 



f(u 



(27) 



The flux-formula makes use of the complete char- 
acteristic information of the system, eigenvalues (char- 
acteristic speeds) and right and left eigenvectors 
through the viscosity term Q. This information is used 
to provide the appropriate amount of numerical dis- 
sipation to obtain accurate representations of discon- 
tinuous solutions without excessive smearing, avoid- 
ing, at the same time, the growth of spurious numeri- 
cal oscillations associated with the Gibbs phenomenon. 
Generic expressions for the characteristic speeds and 
right eigenvectors for the general relativistic hydrody- 
namic equations can be found in Font et al (1999). The 
left eigenvectors have been obtained by Ibanez (1998). 

Notice that the numerical flux is computed at the 
cell interfaces (e.g., i ± 1/2) using information from 
the left and right sides. The state variables, u, must 
be accordingly computed (reconstructed) in advance at 
both sides of a given interface out of the cell-centered 
quantities. The computation of these variables deter- 
mines the spatial order of the numerical algorithm and, 
in turn, controls the local jumps at every interface. If 
these jumps are monotonically reduced the scheme pro- 
vides more accurate initial guesses for the (either exact 
or approximate) solution of the local Riemann prob- 
lems. A wide variety of cell reconstruction procedures 
is available in the literature. As this issue is particularly 
relevant for the simulations we report below, we review 
the most standard choices in the following section. 

In addition, our code has the ability of using dif- 



© RAS, MNRAS 000, 



Nonlinear hydrodynamical evolution of rotating relativistic stars 5 



ferent approximate Riemann solvers: the Roe solver 
(Roe 1981), widely employed in classical fluid dynam- 
ics simulations, with arithmetically averaged states, the 
HLLE solver (Harten, Lax & van Leer 1983; Einfeldt 
1988) and the Marquina solver (Donat & Marquina 
1996), which has been extended to relativistic hydro- 
dynamics by Donat et al (1998). We have performed 
a detailed comparison of the different solvers finding 
good overall agreement among them. Hence, for the fi- 
nal computations reported here we have employed Mar- 
quina's scheme. 

Further specific issues of the hydrodynamical code 
are in order: 1) the reconstructed (left and right) vari- 
ables are the physical (primitive) variables w, except 
for the internal energy density, e, as we use a zero- 
temperature EOS. From these, the remaining (extrap- 
olated) variables are computed algebraically. 2) As we 
are considering adiabatic evolutions, we only solve for 
the first four equations of system ([li]). The internal 
energy (proportional to the rest-mass density) is ob- 
tained algebraically using the EOS. 3) Once all vari- 
ables have been reconstructed, the set of local Riemann 
problems (as a result of the discretization on a numer- 
ical grid) is solved using Marquina's flux-formula. 4) 
The numerical fluxes are computed independently for 
each direction and the time update of the state- vector 
u is done simultaneously using a method of lines in 
combination with a second-order (in time) conserva- 
tive Runge-Kutta scheme, as derived by Shu & Osher 
(1989). 5) After the update of the conserved quanti- 
ties, the primitive variables must be reevaluated. As 
the relation between the two sets is not in closed al- 
gebraic form, the update of the primitive variables is 
done using an iterative Newton- Raphson algorithm. 

4.2 Cell reconstruction procedures 

We turn now to describe the different ways we have 
considered to compute the u L ' R states at every side of 
a cell interface. This is commonly referred to as the 
cell-reconstruction process. In the code we have three 
basic ways of performing such cell reconstruction. 

4.2.1 MUSCL 

First, we can use a second-order TVD scheme (Harten 
1984) of the MUSCL type (Monotonic Upstream 
Schemes for Conservation Laws, van Leer 1979). The 
total variation of a given quantity {u™}"^.^ on a nu- 
merical grid is defined as: 

oo 

7V(u n )= K+i-»?l- (28) 

i= — oo 

A 2-level scheme is called TVD if 

TV(u n+1 ) < TV(u"). (29) 



The TVD property guarantees the suppresion of spu- 
rious oscillations near discontinuities. 

The code uses slope-limiter methods to construct 
second-order TVD schemes by means of monotonic 
piecewise linear reconstructions of the cell-centered 
quantities. In this scheme, uf and \if +1 are computed 
to second-order accuracy as follows: 

uf = u, + Ot{x l+ 1 - Xi) (30) 

uf +1 = u i+ i + a i+ i(x i+ i — x i+ i) (31) 

(x denotes a generic spatial coordinate. In the code it 
can be either r or 8 depending on the direction the 
integration takes place). Our choice for the slope is 
the standard minmod slope which provides the desired 
second-order accuracy for smooth solutions, while still 
satisfying the TVD property: 

a, = mmmod {—^—, ~^^J > ( 32 ) 

where the minmod function of two arguments is defined 
by: 

{a if \a\ < \b\ and ab > 
b if \b\ < \a\ and ab > 
if ab < 

Notice that setting cr* = provides a first-order 
piecewise constant reconstruction scheme. This is then 
equivalent to the original scheme by Godunov (1959). 

4.2.2 ENO 

Our second choice involves the use of high-order (up 
to third-order) ENO schemes (Harten & Osher 1987; 
Harten et al. 1987). In these schemes, given Ui, the cell 
average of a piecewise smooth function, one must con- 
struct a piecewise polynomial function of x of uniform 
polynomial degree r — 1, R(x; u), of the form: 

r-l 

R(x;uj) = y f h,i{x - Xj) 1 (33) 

with x € \xi-i/2, Xj+i/a] and bi t i being the correspond- 
ing polynomial coefficients. This reconstruction poly- 
nomial is essentially non-oscillatory in the sense that it 
satisfies 

TV(R(x-u n+1 )) < TV(R(x-u n )) + 0(Ax r ). (34) 

This property allows the ENO schemes to mantain the 
same order of accuracy at local extrema (whereas the 
TVD schemes always drop to first-order). The left and 
right states are then given by 

uf = R(x t+1/2 ;ui), (35) 

uf +i = iJ(as i+1 /2;ui + x). (36) 

The way the coefficients bi t i are computed for different 
order reconstruction polynomials is explained in the 
appendix of Harten et al (1987). 
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4.2.3 PPM 

Finally, we also use the third-order Piecewise Parabolic 
Method (PPM) of Colella & Woodward (1984). In 
its original design the PPM combined the use of a 
parabolic reconstruction procedure with Godunov's ex- 
act Riemann solver for (Newtonian) ideal gases. Here, 
we use the reconstruction approach in combination 
with Marquina's flux formula. 

In the PPM scheme the interpolated interface val- 
ues Uj + i/2 are obtained using the quartic polynomial 
uniquely determined by the five zone average values 
Ui_2, • • • , Ui+2 in the following way: 



i+1/2 



+ 



= Ui + 



Axi + Axi+i 



(Uj+i - Ui) 



/ 2Ax i+l +Ax t 
{ Ax i + Ax i+1 X(Ul+1 
Axi-i + Axi 

2Ax r + AXi+1 

Axj+i + Ax i+2 
Axi + 2Ax i+ i 



Axi 8 m u i+1 
Axi+i S m Ui f . 



with 
X = 



Axi-i + Axi 



Axi+2 + Ax i+1 



2Axi + Ax t+ i 2Ax t+ i + Axi 



(37) 



(38) 



and 8 m u t = min(|<5ui|, 2|u» - Uj_i|,2|u i+1 - u;|) x 
sign(<5ui) if (u; + i — u»)(ui — u;_i) > and other- 
wise. Additionally 



S\ii 



Axi 



Axi-i + Axi + Axi+i 
2Axi-! + Axi 



+ 



Axi+i + Ax, 

Axi + 2Axi+i 
Axi-i + Axi 



(Uj+1 - Ui) 



(Ui 



(39) 



This algorithm provides a third-order accurate repre- 
sentation of Uj + i/2. For smooth flow solutions the left 
and right values are then given by: 



L 

u i+1 



H+i 



(40) 



In the presence of discontinuities these values are 
modified to ensure the monotone character of the inter- 
polation parabola. Colella and Woodward (1984) pro- 
posed the following modifications: uf = uf = u; if 
(uF-UiXui-uf) < n,uf = 3u<-2uf if(uf-uf)( Ui - 
0.5(uf + uf )) > (uf - uf ) 2 /6 and uf = 3u; - 2uf if 
-« - uf )(ui - 0.5(uf + uf )) > « - uf ) 2 /6. 

Additionally, the PPM scheme incorporates a spe- 
cial treatment of contact discontinuities (contact steep- 
ening) as well as a "flattening" procedure to avoid spu- 
rious post shock oscillations. Details on these more spe- 



cialized aspects of the PPM scheme can be found in the 
original reference of Colella & Woodward (1984). 



4.3 Marquina's flux formula 

The way the numerical flux at a given interface sepa- 
rating the states u L and u R is computed in Marquina's 
solver, is as follows: We first compute the sided local 
characteristic variables and fluxes 



= I p (ul) • u L 
= F(u R ) ■ u R 



= l p (uz,) • f( Ui ), 
= l> fl ).f(Ufl), 



for p = 1, . . . , 5. Here V(ul), F(ur), are the 
(normalized) left eigenvectors of the Jacobian ma- 
trices of the system. Let Ai(uz,), . . . , A 5 (ul) and 
Ai(ur), . . . , As(ur) be their corresponding eigenvalues. 
For k — 1, ... ,5 the procedure is the following: 

• If Afc(u) does not change sign in [ul, ur], then the 
scalar scheme is upwind 



If \ k (u L ) > then 

1 k ik 

4>+ = <Pl, 



= 0, 



else 



4>l = o, 

ik ik 

9- = 9r, 
endif 

• Otherwise, the scalar scheme is switched to the 
more viscous, entropy-satisfying, local-Lax-Friedrichs 
scheme 



a k = max |A fe (u)|, 
uer(ui,,u R ) 



•5(0! + 



ctk^i ), 

k k\ 
R - Uk^r), 



T(ul, ur) is a curve in phase space connecting ul and 
u_r. In addition, a k can be determined as 

a fc = max{jA fc (u L )|, jA fc (u fl )|}. 

Marquina's flux formula is then: 



(41) 



P =i 



where t p (u_l), r p (un), are the right (normalized) eigen- 
vectors of the system. For further details on this solver 
we refer the reader to Donat & Marquina (1996). 



4.4 Numerical grid and boundary conditions 

The code uses spherical polar coordinates (r, 8, <f>) and 
assumes axisymmetry, i.e., there are no <j> derivatives. 
The computational domain in the radial direction ex- 
tends from to 1.2 times the radius of the star. In 
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the polar direction the domain extends from (pole) 
to 7r/2 (equator). For the study of radial pulsations of 
spherical stars we use fine grids of, typically, 400 ra- 
dial grid-points and only one angular grid-point. For 
nonradial pulsations, for which the 9- velocity does not 
vanish, the typical grid-sizes we use are 80 2 to 120 2 . 
For studying quasi-radial pulsations in rotating stars, 
we typically use more points in the radial direction than 
in the angular direction. 

We use a number of boundary grid-zones which 
depends on the stencil size of the different schemes. 
Hence, for the MUSCL and EN02 schemes, we need 
to impose boundary conditions in two additional zones 
at each end of the domain. For the EN03 and PPM 
schemes, the number of additional zones is four. The 
boundary conditions are applied to po,v r ,ve and Vt/> 
and they are as follows: at the center (r = 0) the ra- 
dial fluxes vanish. Hence, the radial velocity is anti- 
symmetric across the origin. The remaining variables 
are symmetric. At the outer radial boundary, the "at- 
mosphere" is re-set to the initial data after each time 
step. At the pole (6 = 0) and equator (6 — n/2), all 
variables are symmetric, except for vg, since the ini- 
tial data have equatorial-plane symmetry and, in this 
paper, we assume that ve is a sum of "polar" (even) 
vector harmonics. 



5 CODE TESTS 
5.1 Shock tube test 

We begin by testing the behaviour of the code and of 
its different schemes in a standard shock tube problem 
in flat spacetime. Although this is a known test and 
the results are, to some extent, well documented in the 
literature (see, e.g., Donat et al 1998, Marti and Muller 
1996) our motivation to include them here is to show 
the correct implementation of the different schemes in 
the present code, by comparing to a non trivial exact 
solution with all types of nonlinear waves: a shock, a 
contact discontinuity and a rarefaction wave. 

Results for the shock tube test are plotted in 
Fig |l|. The initial discontinuous conditions correspond 
to problem 1 in Marti and Muller (1996). The different 
panels show the final state, at t = 0.4, for the veloc- 
ity, density and pressure for all numerical schemes. The 
solid lines represent the exact solution and the symbols 
indicate the numerical solution. We use a grid of 200 
zones spanning a domain of unit length. The density 
and pressure are scaled by a factor of 10 and 20 re- 
spectively. From this figure, we conclude that the im- 
plementation of all schemes is correct as all of them 
reproduce the exact solution and give the correct wave 
speeds. The higher order methods such as EN03 and 
PPM give, as expected, the best results. 




1.2 




-0.2 

0.0 0.2 0.4 0.6 0.8 1.0 

tube length 



Figure 1. Results for the shock tube test for the different 
implemented numerical schemes. The solid lines indicate the 
exact solution, at t = 0.4. The symbols indicate the numeri- 
cal solution for the velocity (circles), density (triangles) and 
pressure (plus signs). The density and pressure are scaled 
by a factor of 10 and 20 respectively. 



5.2 Spherical stars 

Next, we test our code in the non-rotating limit. We 
study the hydrodynamical evolution of non-rotating 
neutron star models in stable equilibrium. Since our 
axisymmetric code uses spherical polar coordinates, it 
can also be run in essentially 1-D, by using only one 
grid point in the angular direction. We determine ra- 
dial pulsation frequencies and compare to frequencies 
expected from perturbation theory. In 2-D, we study 
non-rotating stars that are perturbed from equilibrium 
by an axisymmetric quadrupole perturbation and also 
compare the pulsation frequencies to the ones obtained 
with a perturbative code. 

We will focus our attention on two representative 
neutron star models: Model 1 is a N = 1.5 relativis- 
tic soft polytrope with M/R = 0.056, while model 2 is 
a more relativistic and stiffer polytrope with M/R = 
0.15. Table ^ lists the polytropic index N, polytropic 
constant K, central density p c , mass M and circum- 
ferential radius R for the two models considered (all 
quantities are dimensionless, using c = G = M© = 1) . 
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Initial models of spherical stars 



N K 



p c M R 



Model 1 1.5 4.349 8.10 xlO" 4 0.57 10.11 
Model 2 1.0 100 1.28 xlO" 3 1.40 9.59 



Table 1. Dimensionless (c = G = Mq = 1) equilibrium 
properties of the two representative spherical neutron star 
models. 



5.2.1 1-D evolutions: radial pulsations 

The numerical evolution of initially static, non-rotating 
stars is influenced by the size and sign of the truncation 
errors of the hydrodynamical scheme. We observe the 
following characteristics, when using HRSC schemes: 

(i) The truncation errors at the surface of the star 
initiate small-amplitude radial pulsations. 

(ii) The radial pulsations are dominated by a set of 
discrete frequencies, which correspond to the normal 
modes of pulsation of the star. 

(iii) The numerical viscosity of the finite-difference 
scheme damps the pulsations and the damping is 
stronger for the higher frequency modes. 

(iv) The presence of a constant-density atmosphere 
affects the finite differencing at the surface grid-cells, 
which increases the numerical damping of pulsations 
and also can cause the star to drift to larger densities. 

Thus, we find that any hydrodynamical evolu- 
tion of neutron star models, with the present HRSC 
schemes, will be accompanied by small-amplitude ra- 
dial (or quasi-radial for rotating stars) pulsations. The 
initial amplitude of the radial pulsations and the small 
drift in density converge to zero at the expected or- 
der rate, with increasing resolution. The value of the 
density in the atmosphere region can have a large ef- 
fect on the damping of the pulsations, if it is too large. 
To minimize this effect, we typically set the density of 
the atmosphere equal to 10 -6 times the density of the 
central density of the star. 

Fig. ^ shows the evolution of the radial velocity 
v r at a radius of 0.25-R, with Model 1 as initial model. 
We use the second-order MUSCL scheme with 400 ra- 
dial grid-points. The radial velocity is initially a very 
complex function of time. As we will show next, the 
pulsation consists mainly of a superposition of normal 
modes of oscillation of the fluid. The high frequency 
normal modes are damped quickly and after 20ms the 
star pulsates mostly in its lowest frequency modes. Be- 
cause these oscillations are initiated by the truncation 
errors, the magnitude of the radial velocity is extremely 
small. The velocity oscillates around a non-zero resid- 
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Figure 2. Time evolution of the radial velocity of a spher- 
ical star (model 1). The radial pulsations arc initiated by 
the truncation errors of the finite-difference scheme and are 
mainly a superposition of normal modes of the star. 



ual, which also converges to zero as second order with 
increasing resolution. 

The small- amplitude radial pulsations in the non- 
linear, fixed spacetime evolutions correspond to linear 
normal modes of pulsation in the relativistic Cowling 
approximation (McDermott, Van Horn & Scoll 1983), 
in which perturbations of the spacetime are ignored. 
A Fourier transform of the density or radial velocity 
evolution can be used to compute the normal mode 
frequencies. Fig. ^ shows the Fourier transform of the 
radial-velocity evolution shown in Fig. |^. The normal 
mode frequencies stand out as sharp peaks on a con- 
tinuous background. The width of the peaks increases 
with frequency. The frequencies of radial pulsations 
identified from Fig. |§] are shown in Table |^. 

To compare the obtained frequencies to linear nor- 
mal mode frequencies, we use a code that solves the 
linearized relativistic pulsation equations for the stel- 
lar fluid in the Cowling approximation, as an eigen- 
value problem. In Table |2j we present the results of 
this comparison. The typical agreement between fre- 
quencies computed by the two methods is better than 
0.5% for the fundamental F-mode and the lowest fre- 
quency harmonics Hi — Hi and better than 0.8% for 
the the higher harmonics Hz, — Hg. We note that these 
frequencies are global, as they should be for normal 
mode pulsations, i.e. the frequencies are the same at 
any point in the star, in both the radial velocity and 
density evolution. This is a strong test for the accu- 
racy of the evolution code and our results can be used 
as a testbed computation for other relativistic multi- 
dimensional evolution codes. 

While all four schemes give essentially the same 
radial pulsation frequencies there are some differences 
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Figure 3. Fourier transform of the evolution of the radial 
velocity in Fig. |^| The frequencies are in excellent agree- 
ment with linear normal mode frequencies computed with 
an eigenvalue code (vertical dotted lines). The units of the 
vertical axis are arbitrary. 



Figure 4. Time evolution of the density at half-radius, ini- 
tiated by truncation errors, for a non-rotating initial model 
(model 2). The small secular drift in density, when using 
second order schemes (MUSCL and EN02), is significantly 
reduced or nearly eliminated, using third order schemes 
(EN03 and PPM). 



Radial Pulsation Frequencies. Model 1 



Mode nonlinear (kHz) Cowling (kHz) difference 



F 


1.703 


1.697 


0.3% 


Hi 


2.820 


2.807 


0.5% 


ih 


3.862 


3.868 


0.02% 


Jrh 


4.900 


4.910 


0.2% 


Hi 


5.917 


5.944 


0.4% 


H B 


6.930 


6.973 


0.6% 


H 6 


7.947 


8.001 


0.7% 


H 7 


8.960 


9.029 


0.8% 


H S 


9.973 


10.057 


0.8% 


Hg 


11.030 


11.086 


0.5% 



Table 2. Comparison of small-amplitude radial pulsation 
frequencies obtained with the present nonlinear evolution 
code to linear perturbation mode frequencies in the relativis- 
tic Cowling approximation. The equilibrium model is a non- 
rotating TV = 1.5 relativistic polytrope with M/R = 0.056. 



in the hydrodynamical evolution that are worth to be 
emphasized. Fig. ^shows the density evolution, at R/2, 
for Model 2, using all four different numerical schemes. 
The most striking difference is that the secular drift in 
the density is much smaller for the third-order PPM 
and EN03 schemes than for the second-order MUSCL 
and EN02 schemes. With a resolution of 400 radial 
points, the drift is extremely small: after 10ms the den- 



sity has increased by only 0.3%. With the third-order 
schemes, EN03 and PPM, this drift is considerably 
smaller, being practically unnoticeable in the case of 
PPM. Another difference among the schemes, is that 
each one of them excites the various normal modes of 
pulsation with different amplitude in each mode. In 
EN03, the truncation errors at the surface also ex- 
cite some very high-frequency oscillations, in addition 
to low-order normal modes. These high-frequency os- 
cillations, however, do not represent a problem in the 
identification of normal modes in a Fourier transform. 
Finally, the numerical damping rate of the excited pul- 
sations is similar in the four schemes, with EN02 and 
PPM showing the smallest damping rate. 

5.2.2 2-D evolutions: quadrupole pulsations 

Like radial pulsations, small-amplitude nonradial pul- 
sations can be studied with the present evolution code 
and the obtained frequencies can be compared to per- 
turbation results. We find that the truncation errors 
of the finite difference scheme do not excite nonradial 
oscillations with amplitudes large enough that one can 
identify them accurately in an evolution of unperturbed 
initial configurations. Instead, one has to perturb the 
static initial configuration, using an appropriate eigen- 
function for each nonradial angular index I. For ex- 
ample, to excite the lowest-order quadrupolar (/ = 2) 
pulsations in a spherical star, we perturb the vg ve- 
locity component, using an approximate eigenfunction 
with angular behaviour same as the spherical harmonic 
function Y° an d a simple sin(7rr/i?) radial behaviour. 
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Quadrupole Pulsation Frequencies. Model 2 



Mode nonlinear (kHz) Cowling (kHz) difference 



/ 1.852 1.8843 1.7% 

pi 4.095 4.1099 0.4% 

p 2 6.009 6.0351 0.4% 

p 3 7.858 7.8733 0.2% 

PA 9.683 9.6740 0.1% 



Table 3. Comparison of small-amplitude quadrupole (I = 
2) pulsation frequencies, obtained with the present nonlin- 
ear evolution code, to linear perturbation mode-frequencies 
in the relativistic Cowling approximation. The equilibrium 
model is a nonrotating N = 1.0 relativistic polytrope with 
M/R = 0.15 (Model 2). 



The frequencies of the nonradial modes can then be 
identified in a Fourier transform of the time-evolution 
of the vg velocity component. 

Table ^| shows a similar comparison, as in Table |^, 
for the quadrupole (I = 2) pulsations of model 2. Since 
the nonradial modes have to be computed on a two- 
dimensional grid, we cannot use resolutions as high as 
in the 1-D computations. For a grid-size of 120 x 60 
zones, the difference between frequencies computed by 
the two methods is 1.7% for the fundamental /-mode 
and better than 0.4% for the p-modes p\ — P4. For this 
grid-size, frequencies higher than the p4 mode could not 
be computed accurately, because the grid is to coarse 
to resolve their eigenfunctions (higher harmonic eigen- 
functions have a larger number of nodes in the radial 
direction). 

Fig. M shows the evolution of the vg at a radial dis- 
tance of 0.25i?, for model 2, perturbed with an approx- 
imate quadrupolar egenfunction of small amplitude. In 
this evolution, the PPM scheme was used with an (r, 8)- 
grid of 120 x 60 points. The evolution is mainly a sum 
of the lowest-order quadrupolar pulsation modes of the 
star and allows for the accurate identification of the 
I — 2 normal-mode frequencies. 

5.3 Rotating stars 

We now turn to the evolution of initially stationary, 
uniformly rotating neutron stars. In the rotating case, 
we focus on a rapidly rotating model, with angular 
velocity equal to 92% of the mass-shedding (Kepler) 
limit and equation of state and central density pa- 
rameters same as model 2 in Table In these evo- 
lutions, we observe the same qualitative properties 
as for non-rotating stars and an additional important 
property: the rotation law (angular-momentum distri- 
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Figure 5. Time evolution of the vg velocity component 
at 0.25-R, for a perturbed nonrotating initial model (Model 
2) using the PPM scheme and an (r, 0)-grid of 120 X 60 
points. The evolution is mainly a sum of the lowest-order 
quadrupolar pulsation modes of the star, excited using an 
approximate low-amplitude eigenfunction. 



button), changes near the stellar surface, due to the 
truncation errors of the finite difference scheme. This 
can be attributed to mainly two reasons: (i) the veloc- 
ity component v^ of the fluid has a maximum at the 
surface and the second-order TVD scheme, for exam- 
ple, although high-order accurate in smooth regions of 
the solution, reduces to only first-order at extrema of 
the reconstructed variables. ENO schemes, on the other 
hand, retain the order at local extrema. However, as we 
show below, preserving the angular momentum distri- 
bution near the surface of the star, is still problematic 
with such schemes. This points to a second reason, for 
this behaviour, (ii) The code evolves the relativistic 
momenta (Si) and the velocity components (as well 
as the other "primitive" variables) must be recovered 
through a root-finding procedure, which involves di- 
viding by the density (see Marti & Miiller (1996) for 
details of this procedure). At the surface of the star 
(where the density is very small) this contributes to 
obtaining less than second-order accuracy. 

A typical example of the evolution of a rotating 
star is presented in Fig. which shows the evolution 
of the velocity component v$ obtained with the four 
schemes and with a grid-size of 160 x 40 zones in r and 
respectively. Depicted is the initial equilibrium solu- 
tion (solid line) as a function of the radial coordinate 
(in the equatorial plane) and the final configuration, 
after an evolution time of 5ms, which corresponds to 
4 rotational periods. The figure shows that vj, remains 
close to its initial value in the interior of the star. But, 
near the surface, v<f, decreases with time. It is evident 
from this figure that there are significant differences be- 
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Figure 6. Evolution of the velocity component vs of a ro- 
tating star. The profiles show the radial dependence in the 
equatorial plane (9 = tt/2). Dimensionless units are used in 
both axes. 



tween second-order (top panels) and third-order (lower 
panel) schemes. For the same resolution, the third- 
order schemes result in a more accurate preservation 
of the initial rotation law, both in the interior and 
the surface of the star. We note that, in the interior, 
all schemes give satisfactory results, but for the third- 
order schemes, the difference between the initial and 
final solution is negligible. The more accurate preser- 
vation of the rotation law for the third order schemes is 
emphasized in Fig. 0, which shows the change in v,/, at 
the interior (r = 6, top panel) and at the surface of the 
star (r = 9.8, bottom panel), as a function of time, for 
the four difference schemes. The values depicted corre- 
spond to the equatorial plane of the star (6 — tt/2). At 
the interior, the third-order schemes (EN03 and PPM) 
retain accurately the initial rotation law (PPM shows 
a tiny initial jump to a slightly lowest value), while 
the second-order schemes (EN02 and MUSCL) show a 
very small secular drift (0.8% after 5ms) accompanied 
by oscillations. Correspondingly, at the surface, bot- 
tom panel of Fig. |^, all schemes show a nearly linear 
change in v$ , with the third-order PPM scheme having 
the smallest slope. 

In Fig. we show the momentum component of 
the same rotating star for the four numerical schemes 
as in Fig. []. This is one of the conserved quantities 
u directly evolved with the code. As mentioned pre- 
viously, the recovery of the primitive variables is done 
via a root-finding procedure which involves dividing by 
the density. Both quantities, S$ and p, are very close to 
zero around the surface layer and hence the procedure 
is very sensitive to truncation errors. Whereas the four 
schemes give accurate results for the evolution of S^, 
the computation of v$ near the surface of the star is 
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Figure 7. Evolution of at two different radii of a rotating 
star, for the four different schemes: (a) r = 6 (interior), (b) 
r = 9.8 (surface). The values are measured in the equatorial 
plane (9 = tt/2). In the interior of the star the third-order 
schemes (EN03 and PPM) retain accurately the initial val- 
ues while the second-order schemes (EN02 and MUSCL) 
show a very small secular drift (0.8% after 5ms) accompa- 
nied by oscillations. This drift is due to finite-differencing 
truncation errors. Both second-order schemes behave almost 
identically. Near the surface, r = 9.8 (bottom panel), the 
drift is more pronounced, being smallest for the third-order 
PPM scheme. 



only first-order accurate irrespective of the numerical 
method. 

We plot in Fig. ^| the initial and final configura- 
tions of the azimuthal component of the velocity as a 
function of the angular coordinate, 8, at half the ra- 
dius of the star. The solid line represents the initial 
solution and the dashed line corresponds, again, to a 
final solution after 5ms (four rotational periods). All 
schemes give excellent accurate results, most notably 
EN03 and PPM. 

Quasi-radial pulsations of rotating relativistic 
stars have only been studied in the slow-rotation limit 
(but without the assumption of a fixed spacetime) by 
Hartle & Friedman (1975) (see also Datta et al. 1998). 
We can study such pulsations for rapidly rotating neu- 
tron stars, using appropriate radial eigenfunctions to 
excite them. In Fig. hd we show the time evolution of 
the radial component of the velocity (in the equatorial 
plane) at 0.257?, after an initial small-amplitude radial 
excitation of the stationary star studied in this section. 
The quasi-radial pulsations can be followed accurately 
for many dynamical times. 

We note that the computation of modes in a ro- 
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Figure 8. Evolution of the momentum component of a 
rotating star. The profiles show the radial dependence in the 
equatorial plane (8 = 7r/2). Dimensionless units are used in 
both axes. 
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Figure 9. Evolution of the velocity component of a ro- 
tating star. The profiles show the angular dependence at 
half the radius of the star. 

tating star requires an appropritate excitation, even for 
quasi-radial modes, in order to obtain the correspond- 
ing frequencies with good accuracy in a Fourier trans- 
form of hydrodynamical variables. A detailed study of 
quasi-radial modes in rapidly rotating relativistic stars 
of various equations of state and rotation laws (in the 
Cowling approximation), will appear elsewhere. 



6 DISCUSSION 

We have developed a numerical code which integrates 
the equations of general relativistic hydrodynamics to 
study pulsations of rapidly-rotating relativistic stars 
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Figure 10. Evolution of the radial velocity at 0.25i? for 
a rotating initial model (Model 2) using the PPM scheme 
and a (r, 0)-grid of 160 X 40 points. This figure shows the 
quasi-radial pulsations of the star as a result of a external 
perturbation. 



in a fixed background spacetime. The finite-difference 
code is based on a state-of-the-art approximate Rie- 
mann solver (Donat & Marquina 1996). Our axisym- 
metric relativistic hydrodynamical code is capable of 
accurately evolving rapidly rotating stars for many ro- 
tational periods. We find that, for non-rotating stars, 
small amplitude oscillations have frequencies that agree 
with linear, radial and nonradial, normal mode fre- 
quencies in the Cowling approximation. This study has 
been performed using a representative set of second- 
and third-order TVD and ENO numerical schemes. 

Modern HRSC numerical schemes (as the ones 
used in our code), satisfying the "total variation di- 
minishing" property (Harten 1984), are second-order 
accurate in smooth regions of the flow, but only first- 
order accurate at local extrema. In our rotating star 
evolutions we find that this results in a secular drift of 
the angular momentum distribution near the surface of 
the star. Essentially non-oscillatory schemes (Harten & 
Osher 1987), which retain the order at local extrema, 
were also used to investigate the preservation of the ini- 
tial rotation law near the surface, with results similar 
to the second-order TVD methods. The first-order ac- 
curacy obtained at the surface layer irrespective of the 
method employed pointed out to the ill-posedness of 
the primitive variables recovery procedure as the reason 
of the angular momentum loss. The numerical scheme 
which has the smallest change in the rotation law after 
many rotation periods is the third-order PPM scheme. 

It would be interesting to investigate if the compu- 
tational cost of the present code can be reduced with 
the use of surface-adapted coordinates or fixed-mesh 
refinement and also to analyze whether the change of 
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the rotation law per rotation period will be significantly 
smaller in a frame co-rotating with the star. 

The numerical findings reported in this paper are 
important for the study of small-amplitude and non- 
linear oscillations of rotating neutron stars such as /, 
p and r-modes oscillations. 

Having adopted the PPM third-order scheme as 
our preferred choice for studying hydrodynamical pul- 
sations of rapidly-rotating stars, we plan to investi- 
gate axisymmetric pulsations of rotating proto-neutron 
stars, allowing for various rotation laws and equations 
of state. 
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